acc1_cancer_cells = readRDS("./Data/acc1_cancer_cells_15KnCount_V3.RDS")
all_acc_cancer_cells = readRDS("./Data/acc_cancer_cells_V3.RDS")
acc_all_cells = readRDS("./Data/acc_tpm_nCount_mito_no146_15k_with_ACC1_.RDS")
luminal_pathways = c("CHARAFE_BREAST_CANCER_LUMINAL_VS_BASAL_DN","CHARAFE_BREAST_CANCER_LUMINAL_VS_BASAL_UP","CHARAFE_BREAST_CANCER_LUMINAL_VS_MESENCHYMAL_DN","CHARAFE_BREAST_CANCER_LUMINAL_VS_MESENCHYMAL_UP","HUPER_BREAST_BASAL_VS_LUMINAL_DN","LIM_MAMMARY_LUMINAL_PROGENITOR_UP","SMID_BREAST_CANCER_LUMINAL_B_UP" )
# add luminal pathways
luminal_gs = msigdbr(species = "Homo sapiens") %>%as.data.frame() %>% dplyr::filter(gs_name %in% luminal_pathways)%>% dplyr::distinct(gs_name, gene_symbol) %>% as.data.frame()
suffix = "01_25"
data_to_read = ""
source_from_github(repositoy = "DEG_functions",version = "0.2.24")
source_from_github(repositoy = "HMSC_functions",version = "0.1.12",script_name = "functions.R")
source_from_github(repositoy = "cNMF_functions",version = "0.3.58",script_name = "cnmf_function_Harmony.R")
DimPlot(object = acc1_cancer_cells,pt.size = 2)
FeaturePlot(object = acc1_cancer_cells,features = c("TP63","ACTA2","IL12B","CNN1"))
Warning in FeaturePlot(object = acc1_cancer_cells, features = c("TP63", :
All cells have the same value (0) of IL12B.
patient.ident = all_acc_cancer_cells$patient.ident %>% as.data.frame()
patient.ident[,1] = as.character(patient.ident[,1])
patient.ident[patient.ident[,1] == "ACC1",] = "HMSC"
patient.ident[,1] = as.factor(patient.ident[,1])
all_acc_cancer_cells = AddMetaData(object = all_acc_cancer_cells,metadata = patient.ident,col.name = "patient.ident")
all_acc_cancer_cells = SetIdent(all_acc_cancer_cells, value ="patient.ident")
acc_deg <- FindMarkers(all_acc_cancer_cells, ident.1 = "HMSC",logfc.threshold = 1.5,features = VariableFeatures(all_acc_cancer_cells))
enrichment_analysis(acc_deg,background = VariableFeatures(all_acc_cancer_cells),fdr_Cutoff = 0.01,ident.1 = "HMSC",ident.2 = "ACC",show_by = 1)
hallmark_name = "GO_MITOTIC_CELL_CYCLE"
genesets =getGmt("./Data/h.all.v7.0.symbols.pluscc.gmt")
var_features=all_acc_cancer_cells@assays$RNA@var.features
geneIds= genesets[[hallmark_name]]@geneIds
score <- apply(all_acc_cancer_cells@assays$RNA@scale.data[intersect(geneIds,var_features),],2,mean)
all_acc_cancer_cells=AddMetaData(all_acc_cancer_cells,score,hallmark_name)
#filter:
all_acc_cancer_cells_ccFiltered=all_acc_cancer_cells[,all_acc_cancer_cells@meta.data[[hallmark_name]]< 0.3]
min_threshold = min(all_acc_cancer_cells$GO_MITOTIC_CELL_CYCLE)
max_threshold = max(all_acc_cancer_cells$GO_MITOTIC_CELL_CYCLE)
FeaturePlot(object = all_acc_cancer_cells,features = hallmark_name) + ggtitle("Before cc filtering") & scale_color_gradientn(colours = plasma(n = 10, direction = -1), limits = c(min_threshold, max_threshold))
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
FeaturePlot(object = all_acc_cancer_cells_ccFiltered,features = hallmark_name) + ggtitle("After cc filtering") & scale_color_gradientn(colours = plasma(n = 10, direction = -1), limits = c(min_threshold, max_threshold))
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
patient.ident = all_acc_cancer_cells_ccFiltered$patient.ident %>% as.data.frame()
patient.ident[,1] = as.character(patient.ident[,1])
patient.ident[patient.ident[,1] == "ACC1",] = "HMSC"
patient.ident[,1] = as.factor(patient.ident[,1])
all_acc_cancer_cells_ccFiltered = AddMetaData(object = all_acc_cancer_cells_ccFiltered,metadata = patient.ident,col.name = "patient.ident")
all_acc_cancer_cells_ccFiltered = SetIdent(all_acc_cancer_cells_ccFiltered, value ="patient.ident")
acc_deg <- FindMarkers(all_acc_cancer_cells_ccFiltered, ident.1 = "HMSC",logfc.threshold = 1.5,features = VariableFeatures(all_acc_cancer_cells_ccFiltered))
enrichment_analysis(acc_deg,background = VariableFeatures(all_acc_cancer_cells_ccFiltered),fdr_Cutoff = 0.01,ident.1 = "HMSC",ident.2 = "ACC",show_by = 1)
all_acc_cancer_cells = SetIdent(object = all_acc_cancer_cells,value = "patient.ident") #active snn graph
FeaturePlot(object = all_acc_cancer_cells,features = "MYB",label = T)
new.cluster.ids <- c("cancer", #0
"cancer", #1
"CAF", #2
"cancer", #3
"Endothelial", #4
"cancer", #5
"cancer", #6
"CAF", #7
"CAF", #8
"CAF", #9
"cancer", #10
"CAF", #11
"cancer", #12
"cancer", #13
"cancer", #14
"cancer", #15
"cancer", #16
"WBC", #17
"CAF" #18
)
#rename idents:
acc_all_cells = SetIdent(object = acc_all_cells,value = "RNA_snn_res.1") #active snn graph
names(new.cluster.ids) <- levels(acc_all_cells) #add snn graph levels to new.cluster.ids
acc_all_cells@meta.data[["seurat_clusters"]] = acc_all_cells@meta.data[["RNA_snn_res.1"]]
acc_all_cells = SetIdent(object = acc_all_cells,value = "seurat_clusters")
acc_all_cells <- RenameIdents(acc_all_cells, new.cluster.ids)
# divide "cancer" into patients:
cell_types = acc_all_cells@active.ident %>% as.data.frame()
cell_types[,1]<- as.character(cell_types[,1])
cell_types = cbind(cell_types,acc_all_cells$patient.ident) %>% setnames(old = names(.),
new = c('cell_type','patient'))
cell_types[cell_types$cell_type == "cancer",] = cell_types[cell_types$cell_type == "cancer",2]
# hmsc_rows = (startsWith(x = rownames(cell_types),prefix = "ACC.plate2") | startsWith(x = rownames(cell_types),prefix = "ACC1.")) & cell_types[,1] == "cancer"
# acc_rows = !(startsWith(x = rownames(cell_types),prefix = "ACC.plate2") | startsWith(x = rownames(cell_types),prefix = "ACC1.")) & cell_types[,1] == "cancer"
# cell_types[,1][hmsc_rows] = "HMSC"
# cell_types[,1][acc_rows] = "ACC"
#add to metadata:
cell_types[,2] = NULL
cell_types[cell_types$cell_type == "ACC1",] = "HMSC"
acc_all_cells = AddMetaData(object =acc_all_cells ,metadata = cell_types,col.name = "cell.type")
## {-}
cnv_subtypes = as.data.frame(cutree(infercnv_obj_default@tumor_subclusters[["hc"]][["HMSC"]], k = 2))
names(cnv_subtypes)[1] = "cnv.cluster"
rownames(cnv_subtypes) = rownames(cnv_subtypes) %>% gsub(pattern = "-",replacement = "\\.")
infercnv.observations = data.frame(fread(file = "./Data/inferCNV/infercnv.observations.txt"), row.names=1)
Warning in fread(file = "./Data/inferCNV/infercnv.observations.txt") :
Detected 1332 column names but the data has 1333 columns (i.e. invalid file). Added 1 extra default column name for the first column which is guessed to be row names or an index. Use setnames() afterwards if this guess is not correct, or fix the file write command that created the file to create a valid file.
names_to_keep = colnames(infercnv.observations) %in% (colnames(acc1_cancer_cells) %>% gsub(pattern = "_",replacement = "\\."))
infercnv.observations = infercnv.observations[,names_to_keep]
rownames(cnv_subtypes) = rownames(cnv_subtypes) %>% gsub(pattern = "2\\.",replacement = "2_")
rownames(cnv_subtypes) = rownames(cnv_subtypes) %>% gsub(pattern = "3\\.",replacement = "3_")
acc1_cancer_cells = AddMetaData(object = acc1_cancer_cells,metadata = cnv_subtypes)
DimPlot(acc1_cancer_cells,group.by = "cnv.cluster",pt.size = 2,cols =colors)
original_myo_genes = c( "TP63", "TP73", "CAV1", "CDH3", "KRT5", "KRT14", "ACTA2", "TAGLN", "MYLK", "DKK3")
original_lum_genes = c("KIT", "EHF", "ELF5", "KRT7", "CLDN3", "CLDN4", "CD24", "LGALS3", "LCN2", "SLPI" )
calculate_score(dataset = all_acc_cancer_cells,myo_genes = original_myo_genes,lum_genes = original_lum_genes)
correlation of lum score and myo score: -0.45
correlation of lum score and original lum score: 1
correlation of myo score and original myo score: 1
calculate_score(dataset = acc1_cancer_cells,myo_genes = original_myo_genes,lum_genes = original_lum_genes,lum_threshold = 0,myo_threshold = 0)
correlation of lum score and myo score: 0.06
correlation of lum score and original lum score: 1
correlation of myo score and original myo score: 1
Only HMSC cancer cells:
HPV33_P3 = fread("./Data/HPV33_P3.txt",col.names = c("plate","reads")) %>% as.data.frame()
HPV33_P3.df = HPV33_P3 %>% mutate(
plate = gsub(x =HPV33_P3$plate, replacement = "",pattern = "_.*$")
%>% gsub(pattern = "-P",replacement = ".P")
%>% gsub(pattern = "-",replacement = "_",)
)
HPV33_P3.df = HPV33_P3.df %>% dplyr::filter(HPV33_P3.df$plate %in% colnames(acc1_cancer_cells))
rownames(HPV33_P3.df) <- HPV33_P3.df$plate
HPV33_P3.df$plate = NULL
HPV33_P2 = fread("./Data/HPV33_P2.txt",col.names = c("plate","reads")) %>% as.data.frame()
HPV33_P2.df = HPV33_P2 %>% mutate(
plate = gsub(x =HPV33_P2$plate, replacement = "",pattern = "_.*$")
%>% gsub(pattern = "plate2-",replacement = "plate2_",)
%>% gsub(pattern = "-",replacement = "\\.",)
)
HPV33_P2.df = HPV33_P2.df %>% dplyr::filter(HPV33_P2.df$plate %in% colnames(acc1_cancer_cells))
rownames(HPV33_P2.df) <- HPV33_P2.df$plate
HPV33_P2.df$plate = NULL
HPV33 = rbind(HPV33_P3.df,HPV33_P2.df)
acc1_cancer_cells = AddMetaData(object = acc1_cancer_cells,metadata = HPV33,col.name = "HPV33.reads")
FeaturePlot(acc1_cancer_cells,features = "HPV33.reads",max.cutoff = 40)
data = FetchData(object = acc1_cancer_cells,vars = "HPV33.reads")
print(
data %>%
ggplot(aes( x=HPV33.reads)) +
geom_density()
)
hpv33_positive = HPV33 %>% dplyr::mutate(hpv33_positive = case_when(reads >= 10 ~ "positive",
reads < 10 ~ "negative")
)
hpv33_positive$reads = NULL
acc1_cancer_cells = AddMetaData(object = acc1_cancer_cells,metadata = hpv33_positive)
DimPlot(object = acc1_cancer_cells,group.by = c("hpv33_positive"),pt.size = 2)
Warning in grSoftVersion() :
unable to load shared object '/usr/local/lib/R/modules//R_X11.so':
libXt.so.6: cannot open shared object file: No such file or directory
library(reticulate)
#write expression
acc1_cancer_cells = FindVariableFeatures(object = acc1_cancer_cells,nfeatures = 2000)
vargenes = VariableFeatures(object = acc1_cancer_cells)
hmsc_expression = t(as.matrix(GetAssayData(acc1_cancer_cells,slot='data')))
hmsc_expression = 2**hmsc_expression #convert from log2(tpm+1) to tpm
hmsc_expression = hmsc_expression-1
# hmsc_expression = hmsc_expression[,!colSums(hmsc_expression==0, na.rm=TRUE)==nrow(hmsc_expression)] #delete rows that have all 0
hmsc_expression = hmsc_expression[,vargenes]
write.table(x = hmsc_expression ,file = './Data/cNMF/hmsc_expressionData_2Kvargenes.txt',sep = "\t")
from cnmf import cNMF
name = 'HMSC_cNMF_2Kvargenes'
outdir = './Data/cNMF'
K_range = np.arange(3,10)
cnmf_obj = cNMF(output_dir=outdir, name=name)
counts_fn='./Data/cNMF/hmsc_expressionData_2Kvargenes.txt'
tpm_fn = counts_fn ## This is a weird case where because this dataset is not 3' end umi sequencing, we opted to use the TPM matrix as the input matrix rather than the count matrix
cnmf_obj.prepare(counts_fn=counts_fn, components=K_range, seed=14,tpm_fn=tpm_fn)
cnmf_obj.factorize(worker_i=0, total_workers=1)
cnmf_obj.combine()
cnmf_obj.k_selection_plot()
import pickle
f = open('./Data/cNMF/HMSC_cNMF_2Kvargenes/cnmf_obj.pckl', 'wb')
pickle.dump(cnmf_obj, f)
f.close()
from cnmf import cNMF
import pickle
f = open('./Data/cNMF/HMSC_cNMF_2Kvargenes/cnmf_obj.pckl', 'rb')
cnmf_obj = pickle.load(f)
f.close()
selected_k = 4
density_threshold = 0.1
cnmf_obj.consensus(k=selected_k, density_threshold=density_threshold,show_clustering=True)
/sci/labs/yotamd/lab_share/avishai.wizel/python_envs/miniconda/envs/cnmf_env_6/lib/python3.7/site-packages/scanpy/preprocessing/_simple.py:843: UserWarning: Received a view of an AnnData. Making a copy.
view_to_actual(adata)
usage_norm, gep_scores, gep_tpm, topgenes = cnmf_obj.load_results(K=selected_k, density_threshold=density_threshold)
# calculate usage by expression*genes coefs:
gep_scores = py$gep_scores
all_metagenes= py$usage_norm
plt_list = list()
for (i in 1:ncol(gep_scores)) {
top_genes = gep_scores %>% arrange(desc(gep_scores[i])) #sort by score a
top = head(rownames(top_genes),200) #take top top_genes_num
res = genes_vec_enrichment(genes = top,background = rownames(gep_scores),homer = T,title =
i,silent = T,return_all = T)
plt_list[[i]] = res$plt
}
gridExtra::grid.arrange(grobs = plt_list)
plt_list = list()
for (i in 1:ncol(gep_scores)) {
top_genes = gep_scores %>% arrange(desc(gep_scores[i])) #sort by score a
top = head(rownames(top_genes),10) #take top top_genes_num
message(paste("program ",i,"top genes:"))
print(top)
}
program 1 top genes:
[1] "IGKV5-2" "NDUFB7" "LGALS1" "UBL5" "S100A6" "S100A11" "CUTA" "IFITM3" "JTB" "IL17RB"
program 2 top genes:
[1] "EGR1" "C6orf62" "JUNB" "DNAJA1" "ATF3" "IER2" "ERRFI1" "KLF6" "CDKN1A" "MTHFD2"
program 3 top genes:
[1] "RP1-128M12.3" "RP11-374F3.3" "RP11-403A21.3" "RP11-454L9.2" "AC097500.1" "RP11-463H24.4" "RP11-515O17.2" "RP11-536L3.4"
[9] "PALD1" "AL049758.2"
program 4 top genes:
[1] "LTF" "PIGR" "FMO2" "HSD17B2" "MLPH" "PRR15L" "RF00019.219" "RP11-817J15.2"
[9] "CD14" "CLU"
# Make metagene names
for (i in 1:ncol(all_metagenes)) {
colnames(all_metagenes)[i] = "metagene." %>% paste0(i)
}
#add each metagene to metadata
for (i in 1:ncol(all_metagenes)) {
metage_metadata = all_metagenes %>% select(i)
acc1_cancer_cells = AddMetaData(object = acc1_cancer_cells,metadata = metage_metadata)
}
Note: Using an external vector in selections is ambiguous.
ℹ Use `all_of(i)` instead of `i` to silence this message.
ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
This message is displayed once per session.
FeaturePlot(object = acc1_cancer_cells,features = colnames(all_metagenes))
larger_by = 2
message(paste("larger_by = ", larger_by))
larger_by = 2
acc1_cancer_cells = program_assignment(dataset = acc1_cancer_cells,larger_by = larger_by,program_names = colnames(all_metagenes))
selected_k =4
colors = rainbow(selected_k)
colors = c(colors,"grey")
DimPlot(acc1_cancer_cells,group.by = "program.assignment",pt.size = 2,cols =colors)
show cell cycle program:
hallmark_name = "GO_MITOTIC_CELL_CYCLE"
genesets =GSEABase::getGmt("./Data/h.all.v7.0.symbols.pluscc.gmt")
var_features=acc1_cancer_cells@assays$RNA@var.features
geneIds= genesets[[hallmark_name]]@geneIds
score <- apply(acc1_cancer_cells@assays$RNA@data[intersect(geneIds,var_features),],2,mean)
acc1_cancer_cells=AddMetaData(acc1_cancer_cells,score,hallmark_name)
FeaturePlot(object = acc1_cancer_cells,features = hallmark_name)
cc_vs_program2 = FetchData(object = acc1_cancer_cells,vars = c("metagene.2",hallmark_name))
cor(cc_vs_program2[1],cc_vs_program2[2])
GO_MITOTIC_CELL_CYCLE
metagene.2 0.611673
cnv_vs_hpv = FetchData(object = acc1_cancer_cells,vars = c("cnv.cluster","hpv33_positive"))
test <- fisher.test(table(cnv_vs_hpv))
ggbarstats(
cnv_vs_hpv, cnv.cluster, hpv33_positive,
results.subtitle = FALSE,
subtitle = paste0(
"Fisher's exact test", ", p-value = ",
ifelse(test$p.value < 0.001, "< 0.001", round(test$p.value, 3))
)
)
cnv_vs_cnmf = FetchData(object = acc1_cancer_cells,vars = c("program.assignment","cnv.cluster"))
cnv_vs_cnmf = cnv_vs_cnmf %>% dplyr::filter(program.assignment == "1" |program.assignment == "2" )
test <- fisher.test(table(cnv_vs_cnmf))
ggbarstats(
cnv_vs_cnmf, program.assignment, cnv.cluster,
results.subtitle = FALSE,
subtitle = paste0(
"Fisher's exact test", ", p-value = ",
ifelse(test$p.value < 0.001, "< 0.001", round(test$p.value, 3))
)
)
cnmf_vs_hpv = FetchData(object = acc1_cancer_cells,vars = c("program.assignment","hpv33_positive"))
cnmf_vs_hpv = cnmf_vs_hpv %>% dplyr::filter(program.assignment == "1" |program.assignment == "2" )
test <- fisher.test(table(cnmf_vs_hpv))
ggbarstats(
cnmf_vs_hpv, program.assignment, hpv33_positive,
results.subtitle = FALSE,
subtitle = paste0(
"Fisher's exact test", ", p-value = ",
ifelse(test$p.value < 0.001, "< 0.001", round(test$p.value, 3))
)
)
myb_vs_cnv = FetchData(object = acc1_cancer_cells,vars = c("cnv.cluster","MYB"))
myb_vs_cnv $cnv.cluster = as.character(myb_vs_cnv $cnv.cluster )
ggboxplot(myb_vs_cnv, x = "cnv.cluster", y = "MYB",
palette = "jco",
add = "jitter")+ stat_compare_means(method = "wilcox.test",comparisons = list(c("1","2")))
myb_vs_hpv = FetchData(object = acc1_cancer_cells,vars = c("hpv33_positive","MYB"))
myb_vs_hpv $hpv33_positive = as.character(myb_vs_hpv $hpv33_positive )
ggboxplot(myb_vs_hpv, x = "hpv33_positive", y = "MYB",
palette = "jco",
add = "jitter")+ stat_compare_means(method = "wilcox.test",comparisons = list(c("positive","negative")))+ stat_summary(fun.data = function(x) data.frame(y=15, label = paste("Mean=",round(mean(x),digits = 2))), geom="text") +ylab("log2(MYB)")
hpvReads_vs_myb = FetchData(object = acc1_cancer_cells,vars = c("HPV33.reads","MYB"))
corr = cor(hpvReads_vs_myb$HPV33.reads,hpvReads_vs_myb$MYB)
print("correlation of MYB abd hpv33_reads:" %>% paste(corr %>% round(digits = 2)) )
[1] "correlation of MYB abd hpv33_reads: 0.21"
acc1_cancer_cells_data = acc1_cancer_cells@assays[["RNA"]]@data %>% as.data.frame()
acc1_cancer_cells_data = cor(acc1_cancer_cells_data)
annotation = FetchData(object = acc1_cancer_cells,vars = c("orig.ident"))
colors <- c(seq(-1,1,by=0.01))
my_palette <- c("blue",colorRampPalette(colors = c("blue", "white", "red"))
(n = length(colors)-3), "red")
pht1 = pheatmap(acc1_cor,annotation_col = annotation,fontsize = 6,breaks = colors, color = my_palette,show_colnames = F,show_rownames = F)
acc1_cancer_cells_data = acc1_cancer_cells@assays[["RNA"]]@data %>% as.data.frame()
acc1_cancer_cells_data = cor(acc1_cancer_cells_data)
annotation = FetchData(object = acc1_cancer_cells,vars = c("cnv.cluster"))
colors <- c(seq(-1,1,by=0.01))
my_palette <- c("blue",colorRampPalette(colors = c("blue", "white", "red"))
(n = length(colors)-3), "red")
pht1 = pheatmap(acc1_cor,annotation_col = annotation,fontsize = 6,breaks = colors, color = my_palette,show_colnames = F,show_rownames = F)
cnv_vs_plate = FetchData(object = acc1_cancer_cells,vars = c("cnv.cluster","orig.ident"))
test <- fisher.test(table(cnv_vs_plate))
ggbarstats(
cnv_vs_plate, cnv.cluster, orig.ident,
results.subtitle = FALSE,
subtitle = paste0(
"Fisher's exact test", ", p-value = ",
ifelse(test$p.value < 0.001, "< 0.001", round(test$p.value, 3))
)
)
# creat colours for each group
cnv_vs_plate = FetchData(object = acc1_cancer_cells,vars = c("cnv.cluster","orig.ident"))
rownames(cnv_vs_plate) = rownames(cnv_vs_plate) %>% gsub(pattern = "_",replacement = "\\.")
cnv_vs_plate$cnv.cluster = as.character(cnv_vs_plate$cnv.cluster)
cnv_vs_plate = cnv_vs_plate %>% dplyr::rename(plate = "orig.ident")
annoCol <- list(plate = c(ACC1.P3 = "red",ACC.plate2 = "green"),cnv.cluster =c("1"="green","2" = "red"))
pheatmap(infercnv.observations2,cluster_cols = F,cluster_rows = F, show_rownames = F,show_colnames = F, breaks = breaks,color = colorRampPalette(rev(c("darkred", "white", "darkblue")))(15),annotation_row = cnv_vs_plate,annotation_colors = annoCol)
cnmf_vs_plate = FetchData(object = acc1_cancer_cells,vars = c("program.assignment","orig.ident"))
cnmf_vs_plate= cnmf_vs_plate %>% dplyr::filter(program.assignment == "1" | program.assignment == "2")
test <- fisher.test(table(cnmf_vs_plate))
ggbarstats(
cnmf_vs_plate, program.assignment, orig.ident,
results.subtitle = FALSE,
subtitle = paste0(
"Fisher's exact test", ", p-value = ",
ifelse(test$p.value < 0.001, "< 0.001", round(test$p.value, 3))
)
)
plate_3 = cnmf_vs_plate %>% dplyr::filter((program.assignment == 1 & orig.ident == "ACC1.P3") ) %>% rownames()
plate_2 = cnmf_vs_plate %>% dplyr::filter((program.assignment == 2 & orig.ident == "ACC.plate2")) %>% rownames()
cells = list(ACC1.P3 = plate_3,ACC.plate2 = plate_2)
exceptions_plt = DimPlot(object = acc1_cancer_cells, cells.highlight = cells, cols.highlight = c("green","red"), cols = "gray", order = TRUE,pt.size = 2,sizes.highlight = 2,combine = F)
colors = c("green","red","cyan","purple","grey")
program.assignment_plt = DimPlot(acc1_cancer_cells,group.by = "program.assignment",pt.size = 2,cols = colors,combine = F)
metagene.1_plt <- FeaturePlot(object = acc1_cancer_cells,features = c("metagene.1"),combine = F)
metagene.2_plt = FeaturePlot(object = acc1_cancer_cells,features = c("metagene.2"),combine = F)
lst = list(exceptions = exceptions_plt, program.assignment = program.assignment_plt,metagene.1 = metagene.1_plt,metagene.2 = metagene.2_plt)
for (i in 1:length(lst)) {
cat("### ",names(lst)[i]," \n")
print(
lst[[i]]
)
plot.new()
dev.off()
cat(' \n\n')
}
[[1]]
[[1]]
[[1]]
[[1]]
NA
gep_scores_norm = apply(gep_scores, MARGIN = 2, FUN = min_max_normalize)%>% as.data.frame()
gep_scores_norm = sum2one(gep_scores_norm)
all_metagenes = expression_mult(gep_scores = gep_scores_norm,dataset = all_acc_cancer_cells,top_genes = T,z_score = F)
# Make metagene names
for (i in 1:ncol(all_metagenes)) {
colnames(all_metagenes)[i] = "metagene." %>% paste0(i)
}
#add each metagene to metadata
for (i in 1:ncol(all_metagenes)) {
metage_metadata = all_metagenes %>% select(i)
all_acc_cancer_cells = AddMetaData(object = all_acc_cancer_cells,metadata = metage_metadata)
}
FeaturePlot(object = all_acc_cancer_cells,features = colnames(all_metagenes))
all_acc_cancer_cells = program_assignment(dataset = all_acc_cancer_cells,larger_by = 1.25,program_names = colnames(all_metagenes))
original_myo_genes = c("TP63", "TP73", "CAV1", "CDH3", "KRT5", "KRT14", "ACTA2", "TAGLN", "MYLK", "DKK3")
original_lum_genes = c("KIT", "EHF", "ELF5", "KRT7", "CLDN3", "CLDN4", "CD24", "LGALS3", "LCN2", "SLPI")
orig_myoscore=apply(all_acc_cancer_cells@assays[["RNA"]][original_myo_genes,],2,mean)
orig_lescore=apply(all_acc_cancer_cells@assays[["RNA"]][original_lum_genes,],2,mean)
all_acc_cancer_cells = AddMetaData(object = all_acc_cancer_cells,metadata = orig_lescore-orig_myoscore,col.name = "lum_over_myo")
colors = rainbow(4)
colors = c(colors,"grey")
DimPlot(object = all_acc_cancer_cells,group.by = "program.assignment",cols =colors)
DimPlot(object = all_acc_cancer_cells,group.by = "patient.ident")
FeaturePlot(object = all_acc_cancer_cells,features = "lum_over_myo")
dataset <- iris %>%
group_by(Species) %>%
tidyr::nest() %>%
mutate(
.plot = purrr::map(data, ~ ggplot(.x, aes(x = Sepal.Length, y = Petal.Length)) + geom_point())
) %>%
ungroup()
automagic_tabs(input_data = dataset, panel_name = "Species", .output = ".plot", fig.align='center')
unlink("figure", recursive = TRUE)
dataset$.plot[[1]]
dataset$.plot[[2]]
dataset$.plot[[3]]